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ABSTRACT 

Dynamically significant magnetic fields are routinely observed in molecular clouds, with mass-to-flux ratio 
X = {liZs/GyL/B ~ 1 (here E is the total column density and B is the field strength). It is widely believed that 
“subcritical” clouds with A < 1 cannot collapse, based on virial arguments by Mestel and Spitzer and a linear 
stability analysis by Nakano and Nakamura. Here we confirm, using high resolution numerical models that 
begin with a strongly supersonic velocity dispersion, that this criterion is a fully nonlinear stability condition. 
All the high-resolution models with A < 0.95 form “Spitzer sheets” but collapse no further. All models with 
A > 1.02 collapse to the maximum numerically resolvable density. We also investigate other factors determin¬ 
ing the collapse time for supercritical models. We show that there is a strong stochastic element in the collapse 
time: models that differ only in details of their initial conditions can have collapse times that vary by as much 
as a factor of 3. The collapse time cannot be determined from just the velocity dispersion; it depends also on 
its distribution. Finally, we discuss the astrophysical implications of our results. 

Subject headings: star formation 


1. INTRODUCTION 

Molecular clouds evolve under the influence of self-gravity 
so as to condense part of their mass into dense cores and, ul¬ 
timately, stars. The presence of magnetic fields can prevent 
or delay condensation. The possibility was first studied by 
iMestel & Snitzerl ill 9561) . who noted that the magnetic energy 
and the gravitational energy scale in exactly the same way 
with the radius R of the cloud \/R) if flux freezing obtains. 
They argued that there was therefore a critical mass below 
which a cloud threaded by a particular field strength would be 
unable to collapse. 

A more precise but l ess ge neral argument was advanced by 
iNakano & Nakamural (Il978h . who studied the linear theory 
of a self-gravitating, isothermal, equilibrium sheet of plasma 
threaded by a perpendicular magnetic field. They found that 
magnetic fields stabilize the sheet against gravitational col¬ 
lapse if the mass-to-flux ratio is smaller than 1 IItis/G. 

These results motivate the definition of a dimensionless 
mass-to-flux ratio, 

X=2nVG-, ( 1 ) 

B 

where E is the column density of the sheet, and B is the mag¬ 
netic field strength. The exact coefficient used to define A 
depends somewhat on the geometry of the collapse. Here 
we have chosen the coefficient most relevant to the magnetic 
field geometry adopted in this paper, tending to produce thin 
sheets, in agreement with the expectations for magnetically 
supported clouds. Clouds with A > 1 are termed supercriti¬ 
cal, and clo uds with A < 1 are termed s ubcritical. _ 

Both the iMestel & Snitzerl and the iNakano & Nakamural 
models consider exact equilibria. Molecular clouds are far 
from equilibrium, however, with near-virial, highly super¬ 
sonic velocity dispersion. These internal velocities must arise 
from strong turbulence.* Turbulence might change the stabil- 

* The most plausible alternative to turbulence, some type of weakly dissi¬ 
pative ordered flow, does not emerge naturally in any relevant numerical ex¬ 
periments that we are aware of. The mode-mode coupling is always strong. 
Even circularly polarized Alfven waves, which are exact solutions to the com¬ 


ity properties of the cloud, either by compressing a A < 1 flow 
until it collapses, or by providing turbulent support to a cloud 
with A > 1. 

Many works have suggested tha t turbulence could pro¬ 
vide su pport to star-forming clouds. IChandrasekhar & Fermil 
included turbule nt support in their model for interstel¬ 
lar gaseous structures. IMestel & Snitzerl ill9.56h pointed out 
that turbulence tends to decay, and that turbulence of am¬ 
plitude large enough to support a cloud against self-gravity 
would decay especially quickly, although allowing the pos¬ 
sibility that a strong magnetic field might perhaps allow 
longer lived turbulence. The supersonic linewidths observed 
in molecular clouds were attribut ed to radial motion s insid e 
the cloud instead of turbulen ce bv iGoldreich & KwanI ill 9741) . 
IZuckerman & Palmed lll974h argued that if this interpretation 
were true for all clouds where such fluctuations are observed, 
the star formation rate would be too larg e by at least one or- 


197.51) then suggested that 


der of magnitude. lArons & Max 
the observed velocity fluctuations are due to hydromagnetic 
waves. By the late 1980s, this idea was widely accepted (e.g., 
iShu. Adams. & Lizanolll98^ . In the l ate 1990s, however, a 
succession of numerical etyerirnents llMac Low et al.| |199^ 
[Stone,X)striker, & Gammijll998tlGammie & Ostrikei 11996 ) 
strongly suggested that the damping time of turbulence in 
magnetized molecular clouds is close to the dynamical time. 
If one accepts this, then turbulent pressure can be effective in 
supporting self-gravitating clouds only if it is constantly re¬ 
plenished, in which case the support is perhaps more readily 
identified with the stirring mechanism than with the turbu¬ 
lence itself. 

Other work has tended to emphasize the role of turbulence 
in ini tiating gravitational collapse (e.g., iMac Low & Kless^ 
120041) . Regions with a convergent velocity field will natu¬ 
rally tend to collapse sooner than regions with divergent ve¬ 
locity fields. It seems highly likely that some parts of molec¬ 
ular clouds have strongly convergent velocity fields; is this 
ever enough to overcome the stabilizing effects of the mag- 


pressible equations of motion, suffer from a parametric instab ility with a dy¬ 
namical decay rate ISagdeev & Galeevll969tlGoldsteiill978ft . 
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netic field? Can a subcritical cloud be induced to collapse by 
squeezing, or can a supercritical cloud be prevented from col¬ 
lapsing by the introduction of turbulence? The purpose of this 
paper is to investigate these questions using a simple series of 
numerical experiments. 

The plan of the paper is as follows. In §2 we describe the ex¬ 
perimental design, our numerical methods, and the diffusion 
characteristics of our code (based on the ZEUS algorithm). In 
§3 we describe results, including a “fiducial” run, and the in¬ 
fluence of physical and numerical parameters on the outcome. 
§4 summarizes and discusses astrophysical implications. 

2. DESCRIPTION OE NUMERICAL EXPERIMENTS 

We will consider the simplest possible system that can man¬ 
ifest sub/supercritical behavior; a two-dimensional, periodic 
box containing a magnetized, self-gravitating, isothermal gas. 
Since we are interested in studying the effects of turbulence, 
we will introduce a velocity field in the initial conditions 
with statistical properties similar to those found in interstel¬ 
lar clouds. We will then allow the system to evolve for many 
dynamical times, or until it “collapses.” 

Specifically, we consider a square domain in the y — y plane 
of size Lx L. The z direction points out of this plane; no 
quantity depends on z- The initial fluid density is p, and the 
sound speed, which is constant in space and time, is c^. The 
initial field is B = Bx% where Bx is constant. The strength 
of the field can be characterized by Z = In^/GpL/Bx. We 
set the initial value of {B^ = 0, because otherwise in this z- 
independent geometry, asymptotically there is no collapse. 

The initial velocity field is a Gaussian random field with 
zero d ivergence, constructed as in lOstriker, Gammie, & Stone! 
( Il999ft . The initial velocity field has a power spectrum (v^) 
for 2n/L<k <% {2k jL) and k-= 0. This power spec¬ 
trum is consistent, in 2D, with Larson’s Law ~ Z which 
is equivalent to an energy spectrum ~ k^'^ (in 3D Larson’s 
Law implies k^^). The velocity is normalized so that 

the kinetic energy Ek matches the desired value, and 1/3 of 
the kinetic energy is in motions perpendicular to the x — y 
plane of the simulation. 

2.1. Spitzer sheets and simulation units 

In our experiments we will frequently find that matter flows 
along magnetic field lines to form sheets normal to the field. 
These sh eets are given coherence by the self-gravity of the 
medium. ISt)itzed(ll942h was the first to consider the problem 
of the vertical structure of an infinite, self-gravitating, isother¬ 
mal sheet. Spitzer’s solution turns out to be highly useful in 
understanding the evolution of our simulations. Spitzer found 
that the equilibrium density profile of a sheet of surface den¬ 
sity E is p{z) = (E/2//)sech^(z///) where H = Cs'^/{ kGL). 
The corresponding gravitational potential and field are ps = 
2Ci^log(cosh(z///)) andg = —2;rGEtanh(z///). 

Eor any given Spitzer sheet formed during our simulation 
the surface density parameter E corresponds to pLs, where 
Ls is the extent of the region along the fieldlines that a given 
sheet has collected its mass from. At late times we can assume 
that Ls ~ L: most of the mass originally distributed along the 
fieldlines will be collected into the given sheets. This is true 
for the most massive sheets in the supercritical simulations, 
and it is seen even more clearly in the subcritical simulations, 
where at late times in the simulation a single large scale, sta¬ 
ble sheet incorporates most of the mass of the system. In the 
following we will assume Ls~L and 1, = Lp for the Spitzer 


sheets we are largely interested in — those that have collected 
most mass. 

We nondimensionalize our models by setting L— 1, p = 1, 
and Ci = 1. The simulation time unit is therefore L/cs, the 
sound crossing time. In these units, Newton’s gravitational 
constant equals TTnj^, the Jeans length is 1/nj, the peak den¬ 
sity ps of an equilibrium Spitzer sheet of Ls = L is (7rnj)^/2, 
and its half-thickness H equals {nnj)^^. Notice that, because 
of the periodic boundary conditions the Spitzer sheets are 
slightly distorted, but as long as <C 1, the Spitzer solution 
will be approximately correct. 

Most of the simulations presented in this paper have nj = 
3. This implies that the semithickness of the Spitzer sheet is 
0.011, and if we are to resolve this with at least four grid zones 
we need N > 360, where the resolution of our uniform grid is 
N^. This rather stringent resolution requirement explains why 
we have chosen to study the problem in 2D rather than 3D. 

This Spitzer-sheet model is especially useful when A < 1 
inside a largely ordered magnetic field, able to channel the 
flow into sheets, which later might become unstable and col¬ 
lapse, through accretion, collision, and merger. This sheet 
model, however, is not useful where A ^ 1 and collapse is 
unconstrained by the field. 

Simulation units can be converted to dimensional values 
by assuming for illustrative purposes a typical density nn^ = 
lO^cm^^, and a typical temperature T = lOK. Then the 
sound speed q = 0.19kms^^ fixes the unit of speed, and 
from the Jeans length Lj = Cs{k /Gp)^!^ = 1.9pc we obtain 
the unit of length L = njLj = 5.7pc for our standard value 
nj = 3; a Spitzer sheet would then have a peak density of 
nH 2 = 4.4 X 10^ cm^^, and H = 0.064pc. The unit of mass is 
given by pi? = pn?L? = 1.3 x 10^ M©. The unit of time is 
the sound crossing time L = L/cs « 30Myr. A characteristic 
gravitational contraction time is tg = Lj/c* « lOMyr; free-fall 
collapse times are on the order of 0.3f^ « 3 Myr, depending on 
the geometry of the collapse. 

The sound-crossing time of a Spitzer sheet is ~ H/cs « 
0.3Myr. Erom the Spitzer sheet parameters ps and H it 
is possible to define a characteristic “Spitzer” mass Ms — 
{kH)'^1, = c//(G^E) = in dimensionless units, 

with T. = Lp. The factor of K^ is designed to capture the 
mass inside half a wavelength of the shortest unstable mod e 
of the sheet llLedouxlll951t lElmegreen & Elmegreenlll978h . 
Eor parameters typical of a molecular cloud, Ms = l.bM©, 
which is a suggestive result. This may be compared with 
the thermal Jeans mass Mj = pL? = K^?c ?= 
in dimensionless units, with a typical value of Mj = 
49 (r/lOK)^/^ (nH2/10cm-3)-^/2 

2.2. Ambipolar diffusion 

Our simulation utilizes the ideal MHD equations, which 
have some well-known lim itations. Ambipolar diffu sion is ex¬ 
pected to become relevant JKulsrud & Pearce! 19691) at length- 
scales smaller than the damping length for Alfven waves. 
Lad ~ VAtni, where va = B/ffAnp, and tni = l/{Kni), with 
1.9 X lO^^cm^s^' dDraine. Roberge. & Dalgarnnil983ft . 
The number density of ions n; may depend strongly on en¬ 
vironmental factors, such as th e UV illumination and its at¬ 
tenuation by t he cloud material (ICiolek & Mouschoviasll9^ 
IMcKeell989li : it also depends on chemical properties, such as 
the metal abundance in the gas phase. Eor our fiducial mean 
density p, a representative value could be n, ~ 2 X 10 ^cm 
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largely limited by an assumed metal abundance xm ~ 10 
For the typical peak density of a Spitzer sheet ps in our con¬ 
ditions, n, ~ 6 X 10^‘^cm^^ for cosmic-ray dominated ioniza¬ 
tion, and about ten times larger for regions of the cloud that 
are moderately well UV-illuminated. 

The ambipolar diffusion lengthscale Lad can now be com¬ 
pared with the size L of the computational volume, giving an 
estimate of the scale at which ideal MHD stops being a com¬ 
plete dynamical description of the flow. For our mean density 
profile, we find that this lengthscale is Lad ^ T/35, much 
larger than our typical grid spacing L/512, and comparable to 
our typical sheet thickness 2H — L/45. Flowever, ideal MFID 
is still a good description of the most important portions of 
this study; the regions where mass is collected to form dense 
sheets. For p = ps, the lengthscales are Lad ~ L/700 for a 
UV-dark region, and Lad ~ L/7000 for the more illuminated 
case. The decrease of va with density has contributed to this 
effect, together with the larger 

We keep in mind, however, that densities and ionization 
rates vary widely inside clouds, and so the ambipolar diffu¬ 
sion lengthscales and timescales may vary widely. Turbulent 
conditions inside the flow are also expected to increase the im¬ 
portance of ambipolar diffusion, even at larger lengthscales, 
especially near sharp velocity and magnetic gradients. 

2.3. Stopping criterion 

We must fix some criterion for stopping the numerical in¬ 
tegration; when the density in any zone equals or exceeds the 
largest allowed by the Truelove numerical stability condition 
dTruelove et al.l[T99% the run is terminated and classified as 
having collapsed. The Truelove condition requires that the 
local Jeans length be resolved by some algorithm-dependent 
number Nj of grid zones, typically about 4. This requirement 
sets a maximum resolvable density of px = {N/njNj)^ = 
1820(N/512)^ for nj = 3 on our uniform grid of zones. 
We have found that further integration of Truelove-unstable 
models results in large local fluctuations in the density, which 
can produce “explosions” that corrupt the entire computa¬ 
tional domain. Runs that reach t = 2 without violating the 
Truelove condition are classified as stable. 

We have experimented with other collapse detection 
schemes, because the Truelove criterion has the deficiency 
that it is resolution dependent. In the Tables described below 
we report not only the time tj at which the Truelove condition 
is violated, but also the time fio when 1% of the mass exceeds 
10 times the Spitzer density ps. These times are typically 
close to each other, and both can be considered as measures of 
the onset of gravitational instability. The time t\Q has the ad¬ 
vantage of not depending explicitly on numerical resolution, 
but it can be fooled into producing misleadingly short col¬ 
lapse times by strong density fluctuations, particularly when 
the turbulent kinetic energy is large. 

2.4. Numerical methods and tests 

Our simulations are run on a fixed 2D Cartesian g rid, us¬ 
ing the ZEUS algorithm JStone^Norman|[1992a|™ as im - 
plemented for instance in lOstriker. Gammie. & Storiail1999h . 
ZEUS is a numerical algorithm to evolve ideal (non-resistive, 
non-viscous) non-relativistic MHD flows. It is operator- 
split, representing the fields on a (possibly moving) Eule- 
rian staggered mes h. The magnetic field e volution uses con¬ 
strained transport llEvans & Hawlevll 19881) which guarantees 
that V • B = 0 to machine precision, combined with the method 


of characteristics (iHawlev & Stondfl 995ft . which ensures ac¬ 
curate propagation of Alfven waves. ZEUS is explicit in 
time, and so the timestep Af is limited by the Courant con¬ 
ditions. In our problem, usually the most stringent has been 
Af < Ay/va, where va = Bjy/Anp can take very large values 
in density-depleted regions. A numerical density floor, Pfioor, 
has been set to limit density depletion, preventing Af from be¬ 
coming too small; we have directly tested that this tiny non¬ 
conservation of mass by the code does not alter the simulation 
results regarding collapse in any way. The Poisson equation, 
needed to describe self-gravity, i s solved by Eourier tran sform 
methods, using the EETW code (iFrigo & .Tohnsonl2005h . 

Any Eulerian scheme will cause some diffusion of the mag¬ 
netic held with respect to the mass. It is crucial for our exper¬ 
iment that this nonconservation of A be as small as possible. 
The numerical diffusivity of ZEUS is difficult to estimate be¬ 
cause, unlike a physical resistivity, it is flow dependent. An 
empirical approach is therefore required. 

We have studied conservation of A using two distinct meth¬ 
ods. In the first method we initialize a non-self-gravitating 
box using the same initial data as in our main experiments, 
as described above. We evolve the computation to f = 0.5 
and then damp the velocity held exponentially (with timescale 
Uamp = 0.05) until f = 10. If there were no A diffusion the box 
would return to a uniform density, uniform held state. Diffu¬ 
sion changes A, so the final state consists of a unidirectional 
magnetic held with density and held strength varying only 
perpendicular to the held, from which A can be easily mea¬ 
sured. 

In the second method we initialize a non-self-gravitating 
box using the same initial data as in our main experiments, 
but we evolve the computation only to f = 0.5. We then sam¬ 
ple 80 held lines chosen to lie at equal intervals of the vertical 
component of the vector potential (equivalent to lines equally 

spaced in magnetic flux). We then integrate p/ 
along the held line in the x — y plane,^ using linear interpola¬ 
tion to determine p and B at each position, which immediately 
yields A. 

These two methods give nearly identical results. We there¬ 
fore adopt the second method exclusively, since it can be used 
to probe existing numerical data without any additional, ex¬ 
pensive evolution. 

One possible figure of merit for the diffusion in A is Ox /Aq, 
where Ox is the dispersion in sampled values of A at the final 
instant of the simulation, and Aq is the nominal initial value 
(which we call simply A outside of this subsection). Tables 
[2 and 12 show Oxjhi as a function of resolution and of time 
during a single simulation, respectively. The run shown in 
Table 12 has a resolution of 512^. The key points here are 
that Ox!'^ decreases as resolution increases, and that in ev¬ 
ery case ax/Ao is about 10% or less, which suggests that we 
should be able to measure the critical value of A to similar 
accuracy. 

Evidently Oxj'^ is, converging, but as k, rather than 

the expected This may be because of the existence of 
unresolved regions in the flow where most of the diffusion 
occurs, or it may be the result of irreducible “turbulent” dif¬ 
fusion that is present independent of the magnitude of the ef¬ 
fective numerical diffusion. The numerical diffusion is corre¬ 
lated with the amplitude of turbulence. According to Table|2 

^ This is equivalent to integrating p/ {B\ + along the 3D field- 

line. 
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TABLE 1 

Mass-to-flux diffusion in ZEUS, AS A 

FUNCTION OF NUMERICAL RESOLUTION 


N 

t 

O'a/A) 

Amax/Ao 

^rnin/ ^ 

400 

0.1 

9.0% 

1.26 

0.78 

512 

0.1 

8.6% 

1.31 

0.78 

1024 

0.1 

6.1 % 

1.30 

0.87 

2048 

0.1 

4.4% 

1.10 

0.86 

512 

0.2 

8.3 % 

1.20 

0.79 

2048 

0.2 

4.7% 

1.12 

0.85 


TABLE 2 

MASS-TO-FLUX DIFFUSION IN ZEUS, AS A 
FUNCTION OF TIME 


t 

ox/A) 

^max/ ^ 

Anin/A) 

EK{t) 

0 

0 % 

1 

1 

50 

0.01 

0.2% 

1.01 

0.99 

37 

0.02 

0.9% 

1.02 

0.97 

25 

0.03 

2.1% 

1.05 

0.92 

18 

0.04 

3.6% 

1.10 

0.90 

15 

0.05 

6.4% 

1.18 

0.84 

16 

0.06 

8.0% 

1.32 

0.74 

19 

0.1 

8.6% 

1.31 

0.78 

17 

0.2 

8.3% 

1.20 

0.79 

11 

0.3 

9.6% 

1.20 

0.74 

8.5 

0.34 

11.1% 

1.21 

0.71 

5.5 


much of the diffusion occurs very early in the run, when the 
rms velocity is large. 

Another possible figure of merit is the total variation in A, 
that is, |Amax — Amin|/(2Ao), measuring the possible existence 
of localized diffusion events in addition to the overall diffu- 
sivity of the code. We find in Tables^andlJlthat this quantity 
is typically of the size ~ 2.4(7;^, expected in the mean for the 
half-range of a sample of 80 elements randomly taken from a 
Gaussian distribution. However, the distribution of values of 
X might not always be Gaussian, because localized numerical 
diffusion could be important for some fieldlines. In our tables, 
this may be happening when Amax reaches values as large as 
1.3Ao, even at the relatively high resolution of A = 1024. The 
larger total variation of A observed in those simulations sug¬ 
gests a possible risk of masking the subcritical nature of some 
models.^ However, we also find that the total variation of A 
starts to drop at the even higher resolution of A = 2048; nu¬ 
merical resolution seems apparently able to reduce also this 
more local measure of the diffusivity of the numerical code. 

It is worth noting that ambipolar diffusion in nature is likely 
strong enough to dominate the diffusion measured here. For 
our nominal cloud parameters, we have seen in §2.2 that the 
damping length of Alfven waves is of the same order as the 
expected thickness of an equilibrium sheet. Thus our models 
may misrepresent the situation in nature by tying the fluid too 
closely to the magnetic held. The combination of ambipolar 
diffusion and turbulence (which can drive sharp features for 
the ambipolar diffusion to act on) may be a potent driver of 
variations in mass-to-flux ratio in Galactic molecular clouds. 

3. RESULTS 
3.1. Fiducial run 

^ As a possible example of this effect, in Table one model collapses 
despite having A® = 0.9; also the very low resolution models with N = 128 
that collapse for Aq = 0.744. 


As a guide to the dynamics of our numerical experiments, 
we will first describe a “fiducial” run, whose behavior is in a 
sense typical of the other experiments. This run is supercriti¬ 
cal, with A = 1.5, initial Ek = 50 (equivalent to an rms Mach 
number of 10), nj = 3, and A = 512. The panels ofFigures^ 
and|2lshow how condensation proceeds. At first, small density 
concentrations form due to both ram pressure associated with 
the supersonic velocity fluctuations in the initial conditions 
and fluctuations in the magnetic pressure. These later coa¬ 
lesce into larger clumps, typically oriented perpendicular to 
the magnetic fleld.'^ We have seen that these clumps develop 
into fully stable Spitzer sheets in the simulations of subcriti¬ 
cal clouds; here they can be considered also as approximate 
Spitzer sheets, which later come unstable, as the peak den¬ 
sity of these sheets grows. This density growth takes place 
when the clumps merge or collide, and when matter accretes 
from outside the sheet. The largest density of these clumps 
increases as shown in Fig.|3 slowly at the beginning, but very 
steeply close to the end of the run. The energies, on the other 
hand, vary smoothly in time (Fig. 0}, and are not a good pre¬ 
dictor of the time required for instability. 

At a time fio = 0.325, the fraction of matter denser than 
10 times the nominal Spitzer density includes more than 1% 
of the mass pi%(f) > lOps. Not long afterwards, at a time 
ti =0.341 (soon after the last panel in Fig.[3 the peak density 
of the simulation box exceeds the Truelove limit, forcing an 
end to the run. We conclude that the initial state of this fiducial 
run represents an unstable cloud, able to produce dense cores. 
These tw o times are much larger than the linear e-folding time 
found in iNakan 3 ff^ (« 0.035 for A = 1.5); collecting 
matter into the unstable structures takes a longer time than 
the instability process, linear or nonlinear, and dominates the 
total time necessary to achieve instability in the mildly super¬ 
critical clouds. Using our nominal conversion factors from 
simulation to physical units, fio = 9.8Myr. 

3.2. The nonlinear stability criterion 

To discover how precisely the criticality condition was 
obeyed in the numerical experiments, we considered a series 
of runs with Ek = 50 and Ek = 10 while gradually varying 
A. Table|3lists the collapse times for each of these runs. Ev¬ 
idently the criticality condition is very nearly obeyed in the 
numerical evolutions, and there is no evidence of collapse in¬ 
duced by compression. Indeed, given the diffusion of A mea¬ 
sured in §2, it is remarkable (from a numerical standpoint) that 
we are able to reproduce the condition so accurately. Some 
sense of the “error bars” can be obtained by noticing that the 
A = 1 model with Ek = 50 does collapse, while the A = 1.02 
model with Ek — 10 does not. This suggests that the Nakano 
& Nakamura condition is the true nonlinear stability condi¬ 
tion. 

Tables |4] and 13 show a clear trend to make the simulations 
shorter lived as A increases. This trend is expected from the 
already observed stability criterion.^ Models with A « 1 can 
be quite long lived; the A = 1, Ea" = 50 model persists until 
t = 0.756, or about 23Myr for our nominal cloud parameters. 

The run with A = 1.05 and Ek =10 has been done twice, 
with different values of the numerical density floor; the col- 

'* Clumps do not tend to orient perpendicular to the field in our three di¬ 
mensional models IGammie et alJ2003) . Those runs had a resolution of 256'^, 
however, and the supercritical runs did not have A as close to 1 as the models 
considered here. 

^ The most strongly supercritical models (A > 10) collapse very quickly, 
in around one free-fall time. 
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lapse times are identical up to reasonable precision. This al¬ 
lows us to trust the runs using the larger density floor, which 
are much more convenient because it allows a larger timestep, 
largely controlled by the maximum value of the Alfven speed 
Bj on the grid. From here on, we will not report the val¬ 
ues of this purely numerical parameter in these simulations. 

3.3. Influence of the turbulence energy and distribution 

We have also investigated the effect of the amplitude and 
structure of the initial velocity field. Table0]shows the results 
from a series of runs with nj = 3 and A = 512. The column 
marked “Seed” is the seed used to initiate the random num¬ 
ber generator used to generate the initial velocity field. Runs 
with the same seed but different initial kinetic energies have 
velocity fields that are linearly proportional to each other. 

The first series of runs with seed = 1 show a monotonic 
increase in the lifetime of the cloud with kinetic energy, con¬ 
sistent with results reported elsewhere (iGammie & Ostril^ 
1 1996t [^striker. Gammie. & Stonelll999t) . The effect is weak 
at low energies but more pronounced once Ek > 50. 

An even larger effect is obtained by changing the structure 
of the initial velocity field, i.e. by changing the seed. We find 
that models that differ only in the initial seed can have col¬ 
lapse times that vary by up to a factor of 3. 

To further explore this effect, we performed a series of runs, 
choosing sixty different values of the random seed used to set 
up the shape of the initial velocity distribution. In this se¬ 
ries, we have fixed A = 1.5, Ek = 50, A = 512, and nj=3. 
The results can be seen in Figure]^ showing a wide distri¬ 
bution of collapse times. The total range of this sample goes 
from fx = 0.228 to fj = 0.667, equal to ~ 7 to 20Myr for 
the typical cloud parameters used in §2. The mean time is 
(fx) = 0.368 (^11 Myr); the median is located at fx = 0.355, 
and the peak near fx = 0.3, showing a moderate asymmetry. 
This asymmetry is more pronounced in the tails: fx < 0.2 is 
not observed in the sample; while a few values of fx > 0.5 (at a 
similar distance from the median but in the opposite direction) 
are present in the distribution, corresponding to clouds lasting 
between ^ 15 and 20Myr before collapse, much longer than 
the mean lifetime value. 

This stochastic variation in cloud lifetime doubtless has a 
counterpart in nature. The origin of this variability is clear: 
almost all velocity variations occur at the largest scales, and 
are driven by just a few Fourier modes. If these modes hap¬ 
pen to have the right amplitude and phase then collapse is has¬ 
tened. If they are unfavorable, then collapse can be delayed 
by as much as 12Myr for our nominal cloud parameters. 

3.4. Influence of numerical parameters 

The influence of the density floor Pfloor has already been 
shown in §3.2 to be fully negligible, provided this floor is not 
unreasonably large. 

Numerical resolution, on the other hand, can be quite rel¬ 
evant. Any serious simulation of condensation in a nearly 
critical cloud must be able to resolve the half-thickness H of 
a Spitzer sheet, requiring at the very minimum N > l/H — 
For our fiducial choice of nj=3, this requires a mini¬ 
mum of A > 89, and more reasonably A > 200; any simula¬ 
tion run at smaller resolution would not be exploring the most 
basic physics of mass condensation. However, this require¬ 
ment does not seem to take care of all the effects of numerical 
resolution. 

Table |3 lists simulations where we have varied the num¬ 
ber A of active zones in the grid on each direction. There is 


TABLE 3 

Models with A close to 1 


A 

Ek 

Pfloor 

?10 


1.1 

50 

io-« 

0.304 

0.316 

1.05 

50 

10-^ 

0.387 

0.500 

1.0 

50 

10-^ 

0.497 

0.756 

0.95 

50 

10-^ 

>2 

>2 

0.9 

50 

lO-'^ 

>2 

>2 

1.1 

10 

10-'’ 

0.314 

0.367 

1.05 

10 

lO-® 

0.702 

0.741 

1.05 

10 

lO-'^ 

0.702 

0.741 

1.02 

10 

10-^ 

>2 

>2 

1.0 

10 

lO-'^ 

>2 

>2 

1.0 

10 

lO-" 

>2 

>2 

0.9 

10 

lO-" 

>2 

>2 


Note. — Parameters kept fixed in 
these runs: nj = 3, = 512, random 

seed = 2. 


TABLE4 

Supercritical and 

SUBCRITICAL MODELS 


A 

Ek 

Seed 

Oo 

?T 

1.5 

100 

1 

0.601 

0.615 

1.5 

70 

1 

0.521 

0.533 

1.5 

50 

1 

0.325 

0.341 

1.5 

20 

1 

0.260 

0.270 

1.5 

10 

1 

0.250 

0.257 

1.5 

1 

1 

0.259 

0.269 

1.5 

50 

2 

0.126 

0.211 

1.5 

20 

2 

0.193 

0.200 

1.5 

10 

2 

0.233 

0.244 

1.5 

50 

3 

0.510 

0.537 

1.5 

20 

3 

0.346 

0.360 

1.5 

10 

3 

0.299 

0.310 

1.5 

1 

3 

0.314 

0.324 

1.2 

100 

1 

0.715 

0.744 

1.2 

50 

1 

0.376 

0.400 

1.1 

100 

1 

0.787 

0.910 

1.1 

50 

1 

0.377 

0.534 

1.1 

100 

2 

0.314 

0.372 

1.1 

50 

2 

0.304 

0.376 

1.1 

10 

2 

0.314 

0.367 

1.1 

50 

3 

0.735 

0.746 

0.8 

100 

2 

0.216 

>2 

0.8 

10 

2 

>2 

> 2 

0.8 

10 

3 

0.664 

>2 

0.8 

1 

3 

>2 

> 2 


Note. — Parameters kept fixed 
in these runs: nj — 3, N — 512. 


a clear tendency for the more resolved simulations to delay 
fx. This was in part expected, as the Truelove limit density 
Px = (A/fijAx)^ depends steeply on A. This is not, however, 
the main reason for the observed trend. Peak densities grow 
very quickly in the neighborhood of the condensation time, 
almost nullifying in most cases the influence of the exact mag¬ 
nitude of Px on the value of fx. The quantity fio is expected 
to be less directly dependent on resolution, having a defini¬ 
tion where A does not appear; it still shows some dependence 
on resolution, closely correlated to the dependence shown by 
fx, and probably due to details of the dynamics being more 
revealed at higher resolutions and to reduced numerical diffu¬ 
sion. 

Most of our simulations were performed at A = 512. For 
comparison purposes, a simulation run with the same initial 
conditions as the fiducial run, but with A = 2048, has fx = 
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TABLE 5 

Supercritical AND subcritical 
MODELS RUN AT DIFFERENT NUMERICAL 
RESOLUTIONS 


X 

Ek 

Seed 

N 

tio 


CO 

50 

1 

512 

0.068 

0.077 

oo 

50 

1 

1024 

0.071 

0.083 

1000 

50 

1 

512 

0.068 

0.077 

10 

50 

1 

512 

0.112 

0.124 

1.5 

50 

1 

256 

0.283 

0.283 

1.5 

50 

1 

400 

0.305 

0.312 

1.5 

50 

1 

512 

0.325 

0.341 

1.5 

50 

1 

1024 

0.409 

0.434 

1.5 

50 

1 

1536 

0.429 

0.453 

1.5 

50 

1 

2048 

0.445 

0.481 

1.5 

20 

1 

512 

0.260 

■ 0.270 

1.5 

20 

1 

1024 

0.263 

0.283 

1.5 

10 

1 

512 

0.250 

0.251 

1.5 

10 

1 

1024 

0.253 

0.273 

1.2 

50 

1 

256 

0.319 

0.319 

1.2 

50 

1 

512 

0.376 

0.400 

1.2 

50 

1 

1024 

0.504 

0.523 

1.1 

100 

2 

512 

0.314 

0.372 

1.1 

100 

2 

1024 

0.337 

0.738 

0.9 

100 

2 

512 

0.079 

0.351 

0.9 

100 

2 

576 

0.079 

>2 

0.9 

nw 

1 

512 

> 2 

■ >2 

0.85 

100 

2 

512 

0.075 

>2 

0.8 

100 

2 

512 

0.216 

>2 

0.8 

5U 

1 

256 

> 2 

■ >2 

0.7 

100 

2 

512 

0.435 

>2 

0.6 

50 

1 

256 

>2 

>2 


Note. — Parameter kept fixed in these 
runs: «j = 3. 


TABLE 6 

Models with different 

VALUES OF nj. 



A 

Seed 

N 

?T 

2.5 

1.2 

1 

512 

0.87 


1.2 

2 

256 

0.44 


1.2 

3 

256 

0.83 


1.1 

1 

256 

0.75 


1.1 

2 

256 

0.48 


1.1 

3 

256 

0.82 


1.1 

3 

512 

1.21 


1.1 

4 

256 

0.49 


0.9 

3 

512 

> 2 

2.0 

1.2 

1 

256 

1.77 


1.2 

1 

512 

1.87 


1.2 

2 

256 

1.00 


1.2 

3 

256 

1.50 


1.15 

2 

256 

1.33 


1.15 

3 

256 

>2 


1.1 

1 

256 

> 2 


1.1 

1 

512 

>2 


1.1 

2 

256 

9.77 


1.1 

2 

512 

>2 


1.1 

3 

256 

>2 


1.1 

3 

512 

>2 


1.1 

4 

256 

>2 


0.9 

3 

512 

> 2 


Note. — Parameter kept fixed 
in these runs: Ek = 50. 


0.445 and fio = 0.481, noticeably larger, but not enough to 
change the qualitative conclusions. 

One of the subcritical models with X = 0.9 and = 512 
(Table 0 is remarkable because the run presents an unex¬ 
pected collapse at tj = 0.351. However, the same table shows 


that either changing the numerical seed, or a moderate in¬ 
crease in numerical resolution is enough to suppress this un¬ 
usual behavior; also a small change in X can suppress this 
apparently purely artificial collapse. 

We have also run a few comparison simulations at = 128. 
This resolution is insufficient to represent equilibrium Spitzer 
sheets, and indeed the models collapse even for X = 0.744. 
This is due to the excessively low resolution: it is just enough 
to accommodate one sheet semithickness per grid zone, and, 
through the Truelove stability criterion, it allows only a nar¬ 
row density range, limited by px = {N/rijN-r)^ = 114p = 
2.6ps, insufficient to accommodate an eventual moderately 
large oscillation in the density ps of the equilibrium sheets 
formed in subcritical simulations. 

In another series of tests, we started the simulations from 
an equilibrium sheet, and let it evolve in the presence of a 
very small perturbation. For the values A = 1.5 and A = 2, 
the perturbation grows linearly with e-folding times equal to 
0.036 ±0.001 and 0.027 ±0.001, quite comparable to the val¬ 
ues predicted by th e linear theory, 0.035 and 0.025, known 
from lNakanol ( ll988l) to a precision of 5%. Forvaluesof A < 1, 
the simulations have remained stable up to a time f « 4 during 
a few linear test runs performed at A = 0.9 at various resolu¬ 
tions. 

3.5. Influence of H] 

We have performed a small set of simulations exploring the 
influence of nj on the instability criterion found (Table|^. For 
nj = 3, we had found instability for A > 1.05; for nj = 2.5, 
the numerical requirement is no more stringent than A > 1.1; 
however, for nj = 2.0, the numerical requirement for instabil¬ 
ity becomes A > 1.15. A weak dependence of the numerica l 
criterion on nj had already been predicted bv iNakamil jl988l) 
in the linear regime. Linearly unstable modes have a mini¬ 
mum critical wavelength; if we require that this wavelength 
must fit inside the computational box size L, we find that the 
instability criterion will be approximately® 

A> [l-2/(W)]^\ (2) 

which for small values of nj can be more stringent than the 
infinite disk valu e A > 1, The re sults in Table are consis¬ 
tent with Eq. El iNakand ifT^ presents for the case of a 
finite disk a still more stringent criterion for linear instability, 

A > [l — 4/(7rnj^)] based on the assumption that two crit¬ 
ical wavelengths should fit inside the computational box. Our 
slightly larger unstable range might be related to the geomet¬ 
ric difference between a finite disk and our periodic boundary 
conditions. Spitzer sheets pull in magnetic field lines during 
their formation and contraction (Fig.EJ; this may also allow 
collapse at smaller wavelengths than expected in a purely lin¬ 
ear theory. 

4. CONCLUSIONS 

Our simulations confirm that the single most important el¬ 
ement in determining the long term gravitational stability of 
turbulent magnetized clouds is indeed the mass-to-flux ratio, 
dividing supercritical from subcritical clouds. The relevant 
coeffici ent is that corresponding to a sheet geometry, as de¬ 
rived bv iNakano & Nakamurallll978h . 

® Eollowing INakand 119881) . certain integrals involving generalized Rie- 
mann zeta functions have been replaced by simpler expressions. These ap¬ 
proximations are excellent inside our range of interest «j > 2. 
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Turbulent energy has comparatively little influence on the 
presence or absence of stability, up to Mach numbers ^10. 
Subcritical clouds will develop density concentrations due to 
this turbulence, but under an ideal MHD regime, the conse¬ 
quent increase in magnetic pressure prevents further collapse. 
However, total turbulent energy has some influence on the 
lifetime of supercritical clouds, especially as the Mach num¬ 
ber becomes large enough (of the order of ^ 7 in these simu¬ 
lations). 

More interesting is the fact that turbulence introduces a 
stochastic element. The collapse time cannot be predicted 
with certainty from physical parameters such as the mass and 
field in the cloud, and the typical energy of the turbulence mo¬ 
tions, because the random distributions of velocity and density 
can change the lifetime by some factor, seen to be of the order 
of 3 in one large sample. The resulting distribution of life¬ 
times has an asymmetric tail of unusually long-lived clouds. 
We suggest that the existence of such a tail may introduce a 
bias in the observed samples of star-forming clouds. Most star 
formation will take place in the more frequent, shorter lived 
clouds, while observations of clouds will tend to focus on the 
fewer longer lived ones. 

We have seen that the numerical resolution requirements 
needed to study cloud collapse are very stringent, and we ex¬ 
pect they will be even more stringent in 3D. There is a ne¬ 
cessity of resolving the possible equilibrium structures, such 
as the Spitzer sheets, which we have seen fully formed in the 
subcritical clouds, and partially formed during the run-up to 
instability of the mildly supercritical ones. The thickness of 
these sheets scale with the number nj of Jeans lengths as nj^^. 
Accommodating a large number rij of Jeans lengths inside the 
computational volume will therefore be numerically challeng¬ 
ing. Increasing nj by only a factor of 2 requires increasing 
the space resolution by a factor of 4. Unless adaptive mesh 
refinement (AMR) is used, this requires increasing the simu¬ 
lation runtime by factors on the order of 64 = 4^ in 2D, and 
256 = 4^ in 3D. We anticipate that AMR will be used in many 
of the successful simulations of core formation in the future. 

Numerical stability, through the Truelove condition, sets a 
maximum density that can be accommodated at a given spatial 
resolution. Shocks in strongly turbulent flows have large com¬ 
pression ratios, sometimes requiring increasing resolution in 
order to distinguish a transient density increase due to a shock 
from an authentically unstable accumulation of mass able to 
form a collapsed object. 

We have seen that artificially enforcing numerical density 
floors, even relatively large ones, on the order of times 
the background density, had almost no influence in the evo¬ 
lution of the collapse. This result is again not surprising, 
because wide regions of small density have little influence 
on the dense, self-gravitating regions that undergo collapse. 
Density floors can significantly speed up ideal MHD simula¬ 
tions, whose Courant timestep is often limited by large Alfven 
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speeds Bj in the least dense regions. 

This work is limited due to the periodic boundary con¬ 
ditions. We believe this may have favored the collection 
of clumps into larger clumps until the instability can take 
place. Some simulations occasionally show fast-moving 
clumps flowing past each other, and later merging once one 
of them returns through the other side of the periodic com¬ 
putational volume. The periodic boundary conditions make it 
plausible that sooner or later, most of the mass in a given field¬ 
line will collect into a single clump, which then can undergo 
instability if its mass is even slightly supercritical. In real 
clouds with ordered magnetic fields, clumps inside the same 
fieldline but moving in opposite directions are not expected to 
merge; however, it is improbable this will apply to all of the 
fieldlines and so we expect that the instability will still take 
place in a similar form, albeit with an additional stochastic 
factor in the cloud lifetime. 

Two-dimensionality is also a limitation of this work. It has 
strongly limited the topological possibilities for the fieldlines; 
it is conceivable that the consequent limitations in motion 
have favored the collection of mass i nto massive sheets and 
other structure s. Observations (e.g.. iGoodman et al.l ll99(i 
Crutc hedl2QQ4l). and 3D s imulations and studies (e.g.. lBas'r] 
2000[lGammie et alJ2003h indeed indicate that sheets aligned 
perpendicular to the magnetic field are not always the pre¬ 
ferred possibility for the long term development of clouds. 
More variety of clump shapes is expected in a 3D study. The 
larger variety in motions allowed by a 3D magnetic field is 
expected to enhance the already observed stochastic effects, 
and perhaps might also delay mass collection into potentially 
unstable s tructures. However, even in 3D, the simulations per¬ 
formed by lOstriker, Stone. & Gammi^ ll2001l) suggest that the 
stability criterion will still be dominated by the mass-to-flux 
ratio. 

In some of our models, artificial numerical diffusion has 
turned an initially uniform mass-to-flux ratio A into a non- 
uniform distribution, sometimes with striking effects on the 
numerical stability. While this has a numerical origin, non- 
uniform distributions of mass-to-flux are also expected on as- 
trophysical grounds. For instance, turbulence provides struc¬ 
tures and shocks with small lengthscales and strong magnetic 
gradients, conditions favorable to a localized, efficient am- 
bipolar diffusion, which can redistribute mass and magnetic 
flux independently. Cloud collisions can also merge together 
portions of gas having different masses and magnetic fields. 
We plan to study directly the physical effect of a non-uniform 
mass-to-flux ratio in our future work. 
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Fig. 1. — Colormaps of p in the fiducial run. Snapshots at times t — 0.01, 0.02, 0.04, 0.06, 0.12, 0.18, 0.24, 0.30, and 0.34. The logaiithmic colorscale goes 
from dark blue (saturating on black) to deep red, corresponding to densities going from O.Olps = 0.4441p to lOps = 444. Ip, where ps is the peak density of an 
equilibrium Spitzer sheet for the given parameters. 
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Fig. 3.— Increase of the maximum mass density Pmax with time in the fiducial run. The run finishes when Pmax = 1820, the Truelove value, at a time t — 0.341, 
slightly beyond the plotted region. Some of the transient peaks shown here could have provoked a numerical instability at a resolution smaller than the value 
N — 5\2 adopted for the fiducial run. 
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Fig. 4.— Variations in time of the turbulent kinetic energy Ek, the total magnetic energy Eg, and (minus) the gravitational energy —Eq. 





0.2 


0.4 0.6 

(dimensionless) 

jes of tj. Sixty simulations have been run with the parameters A = 1.5, Ek — 
fx- This plot was constructed by summing 60 Gaussian profiles (with (7 = 0. 
o 20Myr. 




